* GECM w/ multiple predictors DGP

cscript
* Set working directory
cd "~/Dropbox/Whitten_Kagalwala/TS/JOP/Acceptance/KW_Replication/Multiple_Predictors/GECM"

do "GECM-MP-DGP.do"

// Monte Carlo Analyses

matrix mc_summary_stats = J(7, 104, .)

matrix colnames mc_summary_stats = t phi_y ///
phi_adl phi_adl_t1 phi_adl_power b1_adl b1_adl_t1 b1_adl_power g1_adl g1_adl_t1 g1_adl_power b2_adl ///
b2_adl_t1 b2_adl_power g2_adl g2_adl_t1 g2_adl_power lrm_adl_x1 lrm_adl_x1_t1 lrm_adl_x1_power lrm_adl_x2 lrm_adl_x2_t1 lrm_adl_x2_power phi_gecm phi_gecm_t1 ///
phi_gecm_power b1_gecm b1_gecm_t1 b1_gecm_power g1_gecm g1_gecm_t1 g1_gecm_power b2_gecm b2_gecm_t1 b2_gecm_power ///
g2_gecm g2_gecm_t1 g2_gecm_power lrm_gecm_x1 lrm_gecm_x1_t1 lrm_gecm_x1_power lrm_gecm_x2 lrm_gecm_x2_t1 lrm_gecm_x2_power phi1_adl22 phi1_adl22_t1 ///
phi1_adl22_power phi2_adl22 phi2_adl22_t1 b1_adl22 b1_adl22_t1 b1_adl22_power g1_adl22 g1_adl22_t1 g1_adl22_power b2_adl22 ///
b2_adl22_t1 b2_adl22_power g2_adl22 g2_adl22_t1 g2_adl22_power lrm_adl22_x1 lrm_adl22_x1_t1 lrm_adl22_x1_power lrm_adl22_x2 lrm_adl22_x2_t1 ///
lrm_adl22_x2_power b3_adl22 b3_adl22_t1 b4_adl22 b4_adl22_t1 phi1_adl33 phi1_adl33_t1 phi1_adl33_power phi2_adl33 phi2_adl33_t1 ///
phi3_adl33 phi3_adl33_t1 b1_adl33 b1_adl33_t1 b1_adl33_power g1_adl33 g1_adl33_t1 g1_adl33_power b2_adl33 b2_adl33_t1 ///
b2_adl33_power g2_adl33 g2_adl33_t1 g2_adl33_power lrm_adl33_x1 lrm_adl33_x1_t1 lrm_adl33_x1_power lrm_adl33_x2 lrm_adl33_x2_t1 lrm_adl33_x2_power ///
b3_adl33 b3_adl33_t1 b4_adl33 b4_adl33_t1 b5_adl33 b5_adl33_t1 b6_adl33 b6_adl33_t1

* Create a counter to store results in the matrix
local j = 1

foreach t in 25 50 75 100 150 500 1000 {
			parallel sim, expr(phi_adl = r(phi_adl) phi_adl_t1 = r(phi_adl_t1) phi_adl_power = r(phi_adl_power) b1_adl = r(b1_adl) ///
			b1_adl_t1 = r(b1_adl_t1) b1_adl_power = r(b1_adl_power) g1_adl = r(g1_adl) g1_adl_t1 = r(g1_adl_t1) ///
			g1_adl_power = r(g1_adl_power) b2_adl = r(b2_adl) b2_adl_t1 = r(b2_adl_t1) b2_adl_power = r(b2_adl_power) ///
			g2_adl = r(g2_adl) g2_adl_t1 = r(g2_adl_t1) g2_adl_power = r(g2_adl_power) lrm_adl_x1 = r(lrm_adl_x1) ///
			lrm_adl_x1_t1 = r(lrm_adl_x1_t1) lrm_adl_x1_power = r(lrm_adl_x1_power) lrm_adl_x2 = r(lrm_adl_x2) ///
			lrm_adl_x2_t1 = r(lrm_adl_x2_t1) lrm_adl_x2_power = r(lrm_adl_x2_power) phi_gecm = r(phi_gecm) ///
			phi_gecm_t1 = r(phi_gecm_t1) phi_gecm_power = r(phi_gecm_power) b1_gecm = r(b1_gecm) ///
			b1_gecm_t1 = r(b1_gecm_t1) b1_gecm_power = r(b1_gecm_power) g1_gecm = r(g1_gecm) ///
			g1_gecm_t1 = r(g1_gecm_t1) g1_gecm_power = r(g1_gecm_power) b2_gecm = r(b2_gecm) b2_gecm_t1= r(b2_gecm_t1) ///
			b2_gecm_power = r(b2_gecm_power) g2_gecm = r(g2_gecm) g2_gecm_t1 = r(g2_gecm_t1) g2_gecm_power = r(g2_gecm_power) ///
			lrm_gecm_x1 = r(lrm_gecm_x1) lrm_gecm_x1_t1 = r(lrm_gecm_x1_t1) lrm_gecm_x1_power = r(lrm_gecm_x1_power) ///
			lrm_gecm_x2 = r(lrm_gecm_x2) lrm_gecm_x2_t1 = r(lrm_gecm_x2_t1) lrm_gecm_x2_power = r(lrm_gecm_x2_power) ///
			phi1_adl22 = r(phi1_adl22) phi1_adl22_t1 = r(phi1_adl22_t1) phi1_adl22_power = r(phi1_adl22_power) ///
			phi2_adl22 = r(phi2_adl22) phi2_adl22_t1 = r(phi2_adl22_t1) b1_adl22 = r(b1_adl22) ///
			b1_adl22_t1 = r(b1_adl22_t1) b1_adl22_power = r(b1_adl22_power) g1_adl22 = r(g1_adl22) ///
			g1_adl22_t1 = r(g1_adl22_t1) g1_adl22_power = r(g1_adl22_power) b2_adl22 = r(b2_adl22) ///
			b2_adl22_t1 = r(b2_adl22_t1) b2_adl22_power = r(b2_adl22_power) g2_adl22 = r(g2_adl22) ///
			g2_adl22_t1 = r(g2_adl22_t1) g2_adl22_power = r(g2_adl22_power) lrm_adl22_x1 = r(lrm_adl22_x1) ///
			lrm_adl22_x1_t1 = r(lrm_adl22_x1_t1) lrm_adl22_x1_power = r(lrm_adl22_x1_power)  ///
			lrm_adl22_x2 = r(lrm_adl22_x2) lrm_adl22_x2_t1 = r(lrm_adl22_x2_t1) lrm_adl22_x2_power = r(lrm_adl22_x2_power) ///
			b3_adl22 = r(b3_adl22) b3_adl22_t1 = r(b3_adl22_t1) b4_adl22 = r(b4_adl22) b4_adl22_t1 = r(b4_adl22_t1) ///
			phi1_adl33 = r(phi1_adl33) phi1_adl33_t1 = r(phi1_adl33_t1) phi1_adl33_power = r(phi1_adl33_power) ///
			phi2_adl33 = r(phi2_adl33) phi2_adl33_t1 = r(phi2_adl33_t1) phi3_adl33 = r(phi3_adl33) ///
			phi3_adl33_t1 = r(phi3_adl33_t1) b1_adl33 = r(b1_adl33) b1_adl33_t1 = r(b1_adl33_t1) ///
			b1_adl33_power = r(b1_adl33_power) g1_adl33 = r(g1_adl33) g1_adl33_t1 = r(g1_adl33_t1) ///
			g1_adl33_power = r(g1_adl33_power) b2_adl33 = r(b2_adl33) b2_adl33_t1 = r(b2_adl33_t1) ///
			b2_adl33_power = r(b2_adl33_power) g2_adl33 = r(g2_adl33) g2_adl33_t1 = r(g2_adl33_t1) ///
			g2_adl33_power = r(g2_adl33_power) lrm_adl33_x1 = r(lrm_adl33_x1) lrm_adl33_x1_t1 = r(lrm_adl33_x1_t1) ///
			lrm_adl33_x1_power = r(lrm_adl33_x1_power) lrm_adl33_x2 = r(lrm_adl33_x2) lrm_adl33_x2_t1 = r(lrm_adl33_x2_t1) ///
			lrm_adl33_x2_power = r(lrm_adl33_x2_power) b3_adl33 = r(b3_adl33) b3_adl33_t1 = r(b3_adl33_t1) ///
			b4_adl33 = r(b4_adl33) b4_adl33_t1 = r(b4_adl33_t1) b5_adl33 = r(b5_adl33)  ///
			b5_adl33_t1 = r(b5_adl33_t1) b6_adl33 = r(b6_adl33) b6_adl33_t1 = r(b6_adl33_t1)) /// 
			reps(1000) noisily seeds(12344 12344 12344 12344 12344 12344 12344 12344): ecm, obs(`t') rho_e(0.6) 
			  
			matrix mc_summary_stats[`j', 1] = `t'
			matrix mc_summary_stats[`j', 2] = 0.6
		
 
			* ADL (1,1)
			qui su phi_adl
			matrix mc_summary_stats[`j', 3] = r(mean)
			qui su phi_adl_t1
			matrix mc_summary_stats[`j', 4] = r(mean)
			qui su phi_adl_power
			matrix mc_summary_stats[`j', 5] = r(mean)
			qui su b1_adl
			matrix mc_summary_stats[`j', 6] = r(mean)
			qui su b1_adl_t1
			matrix mc_summary_stats[`j', 7] = r(mean)
			qui su b1_adl_power
			matrix mc_summary_stats[`j', 8] = r(mean)
			qui su g1_adl
			matrix mc_summary_stats[`j', 9] = r(mean)
			qui su g1_adl_t1
			matrix mc_summary_stats[`j', 10] = r(mean)
			qui su g1_adl_power
			matrix mc_summary_stats[`j', 11] = r(mean)
			qui su b2_adl
			matrix mc_summary_stats[`j', 12] = r(mean)
			qui su b2_adl_t1
			matrix mc_summary_stats[`j', 13] = r(mean)
			qui su b2_adl_power
			matrix mc_summary_stats[`j', 14] = r(mean)
			qui su g2_adl
			matrix mc_summary_stats[`j', 15] = r(mean)
			qui su g2_adl_t1
			matrix mc_summary_stats[`j', 16] = r(mean)
			qui su g2_adl_power
			matrix mc_summary_stats[`j', 17] = r(mean)
			qui su lrm_adl_x1
			matrix mc_summary_stats[`j', 18] = r(mean)
			qui su lrm_adl_x1_t1
			matrix mc_summary_stats[`j', 19] = r(mean)
			qui su lrm_adl_x1_power
			matrix mc_summary_stats[`j', 20] = r(mean)
			qui su lrm_adl_x2
			matrix mc_summary_stats[`j', 21] = r(mean)
			qui su lrm_adl_x2_t1
			matrix mc_summary_stats[`j', 22] = r(mean)
			qui su lrm_adl_x2_power
			matrix mc_summary_stats[`j', 23] = r(mean)




			* GECM
			qui su phi_gecm
			matrix mc_summary_stats[`j', 24] = r(mean)
			qui su phi_gecm_t1
			matrix mc_summary_stats[`j', 25] = r(mean)
			qui su phi_gecm_power
			matrix mc_summary_stats[`j', 26] = r(mean)
			qui su b1_gecm
			matrix mc_summary_stats[`j', 27] = r(mean)
			qui su b1_gecm_t1
			matrix mc_summary_stats[`j', 28] = r(mean)
			qui su b1_gecm_power
			matrix mc_summary_stats[`j', 29] = r(mean)
			qui su g1_gecm
			matrix mc_summary_stats[`j', 30] = r(mean)
			qui su g1_gecm_t1
			matrix mc_summary_stats[`j', 31] = r(mean)
			qui su g1_gecm_power
			matrix mc_summary_stats[`j', 32] = r(mean)
			qui su b2_gecm
			matrix mc_summary_stats[`j', 33] = r(mean)
			qui su b2_gecm_t1
			matrix mc_summary_stats[`j', 34] = r(mean)
			qui su b2_gecm_power
			matrix mc_summary_stats[`j', 35] = r(mean)
			qui su g2_gecm
			matrix mc_summary_stats[`j', 36] = r(mean)
			qui su g2_gecm_t1
			matrix mc_summary_stats[`j', 37] = r(mean)
			qui su g2_gecm_power
			matrix mc_summary_stats[`j', 38] = r(mean)
			qui su lrm_gecm_x1
			matrix mc_summary_stats[`j', 39] = r(mean)
			qui su lrm_gecm_x1_t1
			matrix mc_summary_stats[`j', 40] = r(mean)
			qui su lrm_gecm_x1_power
			matrix mc_summary_stats[`j', 41] = r(mean)
			qui su lrm_gecm_x2
			matrix mc_summary_stats[`j', 42] = r(mean)
			qui su lrm_gecm_x2_t1
			matrix mc_summary_stats[`j', 43] = r(mean)
			qui su lrm_gecm_x2_power
			matrix mc_summary_stats[`j', 44] = r(mean)


			* ADL (2,2)
			qui su phi1_adl22
			matrix mc_summary_stats[`j', 45] = r(mean)
			qui su phi1_adl22_t1
			matrix mc_summary_stats[`j', 46] = r(mean)
			qui su phi1_adl22_power
			matrix mc_summary_stats[`j', 47] = r(mean)
			qui su phi2_adl22
			matrix mc_summary_stats[`j', 48] = r(mean)
			qui su phi2_adl22_t1
			matrix mc_summary_stats[`j', 49] = r(mean)
			qui su b1_adl22
			matrix mc_summary_stats[`j', 50] = r(mean)
			qui su b1_adl22_t1
			matrix mc_summary_stats[`j', 51] = r(mean)
			qui su b1_adl22_power
			matrix mc_summary_stats[`j', 52] = r(mean)
			qui su g1_adl22
			matrix mc_summary_stats[`j', 53] = r(mean)
			qui su g1_adl22_t1
			matrix mc_summary_stats[`j', 54] = r(mean)
			qui su g1_adl22_power
			matrix mc_summary_stats[`j', 55] = r(mean)
			qui su b2_adl22
			matrix mc_summary_stats[`j', 56] = r(mean)
			qui su b2_adl22_t1
			matrix mc_summary_stats[`j', 57] = r(mean)
			qui su b2_adl22_power
			matrix mc_summary_stats[`j', 58] = r(mean)
			qui su g2_adl22
			matrix mc_summary_stats[`j', 59] = r(mean)
			qui su g2_adl22_t1
			matrix mc_summary_stats[`j', 60] = r(mean)
			qui su g2_adl22_power
			matrix mc_summary_stats[`j', 61] = r(mean)
			qui su lrm_adl22_x1
			matrix mc_summary_stats[`j', 62] = r(mean)
			qui su lrm_adl22_x1_t1
			matrix mc_summary_stats[`j', 63] = r(mean)
			qui su lrm_adl22_x1_power
			matrix mc_summary_stats[`j', 64] = r(mean)
			qui su lrm_adl22_x2
			matrix mc_summary_stats[`j', 65] = r(mean)
			qui su lrm_adl22_x2_t1
			matrix mc_summary_stats[`j', 66] = r(mean)
			qui su lrm_adl22_x2_power
			matrix mc_summary_stats[`j', 67] = r(mean)
			qui su b3_adl22
			matrix mc_summary_stats[`j', 68] = r(mean)
			qui su b3_adl22_t1
			matrix mc_summary_stats[`j', 69] = r(mean)
			qui su b4_adl22
			matrix mc_summary_stats[`j', 70] = r(mean)
			qui su b4_adl22_t1
			matrix mc_summary_stats[`j', 71] = r(mean)
			
			* ADL (3,3)
			qui su phi1_adl33
			matrix mc_summary_stats[`j', 72] = r(mean)
			qui su phi1_adl33_t1
			matrix mc_summary_stats[`j', 73] = r(mean)
			qui su phi1_adl33_power
			matrix mc_summary_stats[`j', 74] = r(mean)
			qui su phi2_adl33
			matrix mc_summary_stats[`j', 75] = r(mean)
			qui su phi2_adl33_t1
			matrix mc_summary_stats[`j', 76] = r(mean)
			qui su phi3_adl33
			matrix mc_summary_stats[`j', 77] = r(mean)
			qui su phi3_adl33_t1
			matrix mc_summary_stats[`j', 78] = r(mean)
			qui su b1_adl33
			matrix mc_summary_stats[`j', 79] = r(mean)
			qui su b1_adl33_t1
			matrix mc_summary_stats[`j', 80] = r(mean)
			qui su b1_adl33_power
			matrix mc_summary_stats[`j', 81] = r(mean)
			qui su g1_adl33
			matrix mc_summary_stats[`j', 82] = r(mean)
			qui su g1_adl33_t1
			matrix mc_summary_stats[`j', 83] = r(mean)
			qui su g1_adl33_power
			matrix mc_summary_stats[`j', 84] = r(mean)
			qui su b2_adl33
			matrix mc_summary_stats[`j', 85] = r(mean)
			qui su b2_adl33_t1
			matrix mc_summary_stats[`j', 86] = r(mean)
			qui su b2_adl33_power
			matrix mc_summary_stats[`j', 87] = r(mean)
			qui su g2_adl33
			matrix mc_summary_stats[`j', 88] = r(mean)
			qui su g2_adl33_t1
			matrix mc_summary_stats[`j', 89] = r(mean)
			qui su g2_adl33_power
			matrix mc_summary_stats[`j', 90] = r(mean)
			qui su lrm_adl33_x1
			matrix mc_summary_stats[`j', 91] = r(mean)
			qui su lrm_adl33_x1_t1
			matrix mc_summary_stats[`j', 92] = r(mean)
			qui su lrm_adl33_x1_power
			matrix mc_summary_stats[`j', 93] = r(mean)
			qui su lrm_adl33_x2
			matrix mc_summary_stats[`j', 94] = r(mean)
			qui su lrm_adl33_x2_t1
			matrix mc_summary_stats[`j', 95] = r(mean)
			qui su lrm_adl33_x2_power
			matrix mc_summary_stats[`j', 96] = r(mean)
			qui su b3_adl33
			matrix mc_summary_stats[`j', 97] = r(mean)
			qui su b3_adl33_t1
			matrix mc_summary_stats[`j', 98] = r(mean)
			qui su b4_adl33
			matrix mc_summary_stats[`j', 99] = r(mean)
			qui su b4_adl33_t1
			matrix mc_summary_stats[`j', 100] = r(mean)
			qui su b5_adl33
			matrix mc_summary_stats[`j', 101] = r(mean)
			qui su b5_adl33_t1
			matrix mc_summary_stats[`j', 102] = r(mean)
			qui su b6_adl33
			matrix mc_summary_stats[`j', 103] = r(mean)
			qui su b6_adl33_t1
			matrix mc_summary_stats[`j', 104] = r(mean)
			
			* Iterate the counter
			local j = `j' + 1
		
	}


clear
set obs 35
svmat double mc_summary_stats, names(col)
save gecm.dta, replace



save gecm.dta, replace

* Bias calculations
gen bias_adl_phi = phi_adl - phi_y
gen bias_adl_b1 = b1_adl - 2
gen bias_adl_g1 = g1_adl - 1.5
gen bias_adl_b2 = b2_adl - 1
gen bias_adl_g2 = g2_adl - 1.25
gen bias_adl_lrm_x1 = lrm_adl_x1 - ((2 + 1.5)/(1-phi_y))
gen bias_adl_lrm_x2 = lrm_adl_x2 - ((1+1.25)/(1-phi_y))

gen bias_gecm_phi = phi_gecm - (phi_y - 1)
gen bias_gecm_b1 = b1_gecm - 2
gen bias_gecm_g1 = g1_gecm - (2+1.5)
gen bias_gecm_b2 = b2_gecm - 1
gen bias_gecm_g2 = g2_gecm - (1 + 1.25)
gen bias_gecm_lrm_x1 = lrm_gecm_x1 - ((2 + 1.5)/(1 - phi_y))
gen bias_gecm_lrm_x2 = lrm_gecm_x1 - ((1 + 1.25)/(1 - phi_y))

gen bias_adl22_phi = phi1_adl22 - phi_y
gen bias_adl22_b1 = b1_adl22 - 2
gen bias_adl22_g1 = g1_adl22 - 1.5
gen bias_adl22_b2 = b2_adl22 - 1
gen bias_adl22_g2 = g2_adl22 - 1.25
gen bias_adl22_lrm_x1= lrm_adl22_x1 - ((2 + 1.5)/(1-phi_y))
gen bias_adl22_lrm_x2= lrm_adl22_x2 - ((1 + 1.25)/(1-phi_y))
gen bias_adl22_phi2 = phi2_adl22 - 0
gen bias_adl22_b3 = b3_adl22 - 0
gen bias_adl22_b4 = b4_adl22 - 0

gen bias_adl33_phi = phi1_adl33 - phi_y
gen bias_adl33_b1 = b1_adl33 - 2
gen bias_adl33_g1 = g1_adl33 - 1.5
gen bias_adl33_b2 = b2_adl33 - 1
gen bias_adl33_g2 = g2_adl33 - 1.25
gen bias_adl33_lrm_x1 = lrm_adl33_x1 - ((2 + 1.5)/(1-phi_y))
gen bias_adl33_lrm_x2 = lrm_adl33_x2 - ((1 + 1.25)/(1-phi_y))
gen bias_adl33_phi2 = phi2_adl33 - 0
gen bias_adl33_phi3 = phi3_adl33 - 0
gen bias_adl33_b3 = b3_adl33 - 0
gen bias_adl33_b4 = b4_adl33 - 0
gen bias_adl33_b5 = b5_adl33 - 0
gen bias_adl33_b6 = b6_adl33 - 0

save gecm_mp.dta, replace
